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1. Summary 


1.1 Regular Perturbation Analysis 

During this reporting period research was directed at evaluating the regular perturbation method 
described in details in [1]. Closed-loop simulations were performed with a first order correction 
including all of the atmospheric terms. In addition, a method was developed for independently 
checking the accuracy of the analysis and the rather extensive programming required to implement 
the complete first order correction with all of the aerodynamic effects included. This amounted to 
developing an equivalent Hamiltonian for the first order analysis and evaluating it by quadrature. 
The result was compared to the Hamiltonian computed from the first order analysis. A second 
order correction was also completed for the neglected spherical Earth and back-pressure effects. 
Finally, an analysis was begun on a method for dealing with control inequality constraints. 

To date, the results on including higher order corrections do show some improvement for this 
application, however we do not know at this stage if significant improvement will result when the 
aerodynamic forces are included. If the result is negative, then our recommendation is that the 
method of Matched Asymptotic Expansions (MAE) be explored as the next major step in this 
research effort. The results from a parallel research effort on aeroassisted orbit transfer trajectories 
indicate that the regular perturbation analysis under current investigation actually plays the role of 
the inner expansion in a MAE analysis. The outer solution in a MAE analysis provides a correction 
currently not available from a regular expansion. We would like to explore if a similar situation 
holds for the dynamics associated with launch vehicle trajectories. 

1.2 Finite Element Analysis 

The weak formulation for solving optimal control problems has now been extended in order to 
account for state inequality constraints. The formulation has been tested on three example problems 
and numerical results have been compared to the exact solutions. Development of a general- 
purpose computational environment for the solution of a large class of optimal control problems is 
now well underway. An example, along with the necessary input and the output, is given. 

2. Research Accomplishments 

2.1 Regular Perturbation Analysis 
Closed-loop Simulation 

Figures 1 and 2 compare the performance of the closed-loop control solutions generated by two 
different methods, with the open loop optimal solution generated using a multiple shooting 
method. A first order correction was made in each case for the neglected spherical Earth & engine 
back-pressure effects. The simulation results are for lift and drag set to zero. In Method 1, the 
control update interval was 1 second and within each interval the control was held constant. The 
control was determined by repeatedly calculating a new zero order solution and performing a 
quadrature at every update. Method 2 was based on a pre-calculated quadrature for a fixed zero 
order solution corresponding to the conditions near launch. See [1] for details on these two 
Methods. Nearly continuous control updating was used for Method 2 because the computational 
effort is trivial. It amounts to solving a set of 4 linear equations to generate the on-line control. A 
mid-point extrapolation scheme (accuracy equivalent to a Runge-Kutta 7/8) was used in both 
methods for the simulation. Table 1 gives a comparison of the terminal conditions and the 
performance index. 
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Table IrTerminal Values Comparison 


hr 

V f 

Yf 

tf 


Method 1 
148160m 
7857.58m/s 
0.001 deg 
355.612s 


Method 2 

148147m 

7864.99m/s 

0.035deg 

355.744s 


Optimal 
148160m 
7858. 2m/s 
0 

355.591s 


The results show a dramatic improvement in comparison to the open loop solutions reported in [1] 
for these two methods. In [1] the trajectories were obtained from a single calculation at launch, 
and the trajectories were constructed by simply summing the zero order solution and integrated first 
order dynamics. For the results shown here, the solutions were obtained by integration of the 
complete dynamics, with the control computed from the perturbation analysis. 

Figures 3 and 4 show the closed loop simulation results for Method 1 including the aerodynamic 
forces in the first order correction. Since the zero order solution gives an unrealistically high angle- 
of-attack (approximately -45 deg.) at launch, the simulation was started at an altitude of 10525 
meters, so that the zero order solution for alpha was still within the range of the tabulated 
aerodynamic data. Figure 3 clearly indicates the onset of an instability in alpha at this altitude. The 
slight increase in alpha near the end is due to a numerical problem that can be removed at a later 
date. 


It was not known at this point if this instability was due to an analysis and/or programming error, 
or due to the inability of the regular perturbation analysis to account for aerodynamic effects using 
a first order correction. It could also be that a second-order correction would not significantly 
improve matters, since at best we are forming an asymptotic series solution to the problem. Thus 
we decided to develop an independent check on our results before proceeding to a second order 
analysis, which is described in the next section 

Checking the First-order Analysis 

Checking was performed by monitoring a Hamiltonian function which corresponds to the first 
order necessary conditions when viewed as being derived from an equivalent optimization 
problem. This Hamiltonian is different from the first order expansion of the original Hamiltonian 
for the full nonlinear dynamics. The first order Taylor's series expansion of the original 
Hamiltonian does not correspond to the costate equations and the optimality condition of the first- 
order dynamics [2], So a new Hamiltonian ( H ) was derived which has the following form: 

H = + fj°u, + ^ 41 ° + (t - + g o)TX, 

MMio MM 

dx 3u 

3(fU) .To a(flX) 

3u 3u 



Since the first order system is time-varying, the Hamiltonian is not constant. The first order 
analysis is checked by realizing the following two expressions: 
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( 2 . 2 ) 


d// _ d//(x, X , u*) 
dt at 


// = X T f(x, X, u*) (2.3) 

where for the right hand side of (2.2) we mean the partial derivative of the expression in (2.1). 
The Hamiltonian was computed in two ways. First by numerical integration of (2.2) along the 
trajectory with an arbitrary initial condition. Second, by direct substitution of the state, costate and 
control values from the first order solution into (2.3). If the analysis and programming is correct, 
the difference between the two ways of computing H should be stagewise constant. This was 
verified by the results in Figures 5 and 6. In this setting, both the zero and first-order optimality 
conditions and their costate equations were verified because H also depends on the zero order 
solution. The difference in the two calculations is zero to within 4 significant digits. 

Second-or der Correction 

A second order analysis was also carried out to determine if any improvement results in 
comparison to the first order solution. At this stage, the second order analysis including the 
aerodynamic forces is not completed. However, a second order correction for the spherical Earth 
and back-pressure effects was evaluated and the results are depicted in Figures 7-11. These are the 
open-loop histories obtained by summing the forward integration results for each corrected term. 
In integrating the second order dynamics, the first order state and costate histories are required. 
This is done by a forward integration of the first order dynamics using the known initial values for 
Xi(to) and the calculated initial costate correction X](to). The histories are stored for a sufficient 
number of sample points, and retrieved using a piecewise linear interpolation for the integration of 
the second order dynamics. 

In examining the results of Figures 7-11 it should be noted that the pattern throughout is that the 
first-order correction over-corrects the zeroth-order solution, and that the secord-order correction 
over- corrects the first-order solution. Unfortunately, the error is not significantly decreased by the 
second-order correction in most of the results, with the exception of the Xu profile which shows a 
dramatic improvement. The estimates for the initial values of the costates and the final time 
(performance index) are compared in Table 2. 

Table 2. Performance Comparison Of The Open-loop Results 



X v (0)/s 2 m- 1 

A^(0)/s2m-i 

^(OVsm-i 

tf/s 

Zeroth-order 

0.20156e-l 

0.19334e-l 

0.54560e-4 

360.047 

First-order 

0.39188e-l 

0.22036e-l 

0.56468e-3 

354.335 

Second-order 

0.35344e-l 

0.20899e-l 

0.57143e-3 

355.254 

Optimal 

0.37352e-l 

0.20868e-l 

0.60304e-3 

355.606 


Control Inequality Constraint 

Preliminary work on addressing control inequality constraints (C(x,u) < 0) is under investigation. 
This approach makes use of a slack variable (cr) to transform to a strict equality constraint [3]. The 
necessary conditions are as follows: 
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* = f + eg ; t e [to,tf] (2.4) 

X = -flX - C x p - egjX (2.5) 

0 = fj^. + C u p + eg JX (2.6) 

0 = an (2.7) 

0 = C + i<x 2 (2.8) 

Equations (2.4-2.6) are derived from the necessary conditions on the augmented Hamiltonian, 

H = f T A, + eg T A, + (C + ;kx 2 )p (2.9) 


When the trajectory is on the constraint, a = 0. When it is off the constraint, p = 0. Note that the 
product is zero at every instant. Alternatively, (2.7) can be derived if we realize that the slack 
variable can be treated as a control variable and then use the optimality condition H a = 0. 
Equations (2.7) and (2.8) provides the additional information needed to determine a and p. 

To obtain the zeroth- and higher-order formulations, we simply need to carry out the expansions 
including: 


p = po + epi + £?P2 + ... (2.10) 

a = ao + eai + e ?(*2 + ... (2.11) 

Substituting these expansions and equating like powers in e, the algebraic equations (2.6-2.8) can 
be grouped as (note that for simplicity, we consider a scalar u case): 


‘(f u T X)2 + C2 u Po c° o' 

u i 


(f u^)h + CunMd 


r f ° T i 

•u 

0 a 0 Po 

C° 0 a 0 

aj 


0 

c° 

'“K 

»j- 

0 

0 


^j— 1 * ^j— 1 ) 

0 

L 0 

where j = 1, 2, ... 

Solving the control constraint problem requires a guess of the switching structure. This is true of 
all indirect methods. In this case, it is the switching structure on the zero order solution that 
matters. The method requires that the zero order solution captures the true switching structure 
because it affects the matrix on the left-hand-side of (2.12). It is this matrix which subsequently 
produces the control correction that leads to a better approximation. If the matrix is singular the 
method and the expansion technique will fail. On case that does lead to a singular matrix is the 
touch-point switching structure, where 

a 0 = po = 0 


( 2 . 12 ) 
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det< 


(fftfl+Oo CS 0 

0 oc 0 p 0 

C° 0 a 0 

For some simple cases, it may be possible to incorporate the control constraint in the zero order 
problem, thus capturing the true switching structure. However, for the launch vehicle problem, 
where the dynamics are nonlinear and time varying, and incorporation of any form of control 
constraint will make the derivation of an analytic solution difficult. Further analysis is required to 
see whether any simplification is possible. 

2.2 Finite Element Analysis 

2.2.1 Extension of the Analysis. The method based on the weak Hamiltonian formulation 
derived in [4] and [5] has now been extended to handle problems with state inequality constraints. 
An outline of the derivation and a simple example problem are given in Appendix A. (Even more 
details of the derivation can be found in [6], a copy of which will be sent to the Technical Monitor 
as soon as it is complete.) 

The derivation proceeds in the following manner. It is desired to develop a solution strategy for 
optimal control problems with state inequality constraints based on finite elements in time. In an 
attempt to make the solution scheme as general as possible, all strong boundary conditions are 
transformed into natural boundary conditions. This is done so that the shape functions can be 
chosen from a less restrictive class of functions, which enables one to choose the same shape 
functions for every optimal control problem. 

The idea of transforming the strong boundary conditions to natural boundary conditions [7] 
revolves around adjoining a constraint equation to the performance index with an unknown 
Lagrange multiplier. The variation of the performance index is then taken in a straightforward 
manner. Through appropriate integration by parts, it is possible to show that the Euler-Lagrange 
equations are identical to those derived in classical textbooks [8] and that the boundary conditions 
are the same, only stated weakly instead of strongly. 

2.2.2 Development of a General Code. The weak formulation is capable of solving optimal 
control problems that have continuous states, costates, and controls, and problems with 
discontinuities arising from staging (i.e., discontinuities in the system equations), control 
inequality constraints and state inequality constraints. The algebraic equations which come from the 
weak formulation may be derived prior to specifying the problem to be solved. It is this feature in 
particular that allows for a general problem-solving environment to be created. 

The main goal of the general code is to reliably solve a large class of optimal control problems with 
a minimum of user interaction. Specifically, it is desired to create an environment where the user 
does not have to write subroutines. To this end, the general code is being developed on a SUN 
3/260 workstation and requires a FORTRAN 77 compiler, MACSYMA [9], and the Harwell 
subroutine library [10]. The general procedure can be broken into three parts that must interface 
together. The first part is the FORTRAN code. This code contains all the subroutines necessary to 
solve any of the optimal control problems described above. However, if certain problems require 
table look-up routines (such as aerodynamic data for a rocket model), then these subroutines must 
be given by the user and interfaced to the rest of the general code. Thus, there may be a need for 
some user programming for certain problems. The second pan of the general procedure is the use 
of MACSYMA. The user must supply an input file specifying the problem. This input file is in 
symbolic form and will be loaded into MACSYMA. MACSYMA will then evaluate all the 
necessary expressions and automatically generate the FORTRAN code. This code is spliced into a 


) = (fu ^)u a o + Cu 2 p 0 = 0 


(2.13) 
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template file and becomes one of the subroutines. The third and final part of the general procedure 
will consists of subroutines to generate initial guesses that will reliably converge. Homotopy 
methods are the prime candidates for this. A very simple type of homotopy method described in 
[11] is being used. This method converts the algebraic equations to initial-value ordinary 
differential equations. A second-order Runge-Kutta method is used to integrate the equations and 
obtain initial guesses for a Newton-Raphson method. This method has worked on all the problems 
tested to date. 

The general code is still being developed at this time. Currently, the code can handle problems 
with continuous states, costates, and controls, problems with control inequality constraints, and 
problems with state inequality constraints that only touch the constraint boundary. The general 
code is now functional (but not complete) for a large class of optimal control problems. An 
example problem demonstrating the use of the code is given in Appendix B. 

3. Future Research 

In the perturbation analysis area we plan to complete the second order analysis, and to perform 
both open loop and closed loop comparisons to the first order results and to the optimal solution. 
However, we are skeptical at this point that second order correction will remove the instability 
observed in the first order results when aerodynamic forces are included. Along this line we plan 
to spend some time investigating the potential that Matched Asymptotic Expansions has for 
improving the solutions that we have obtained to date. We will also continue investigating the 
control inequality constraint formulation. Results will first be developed for several simpler 
problems to evaluate its potential for application to launch vehicle guidance problems. 

In the finite element analysis area we plan to complete the development of the general code so as to 
make it applicable to all types of optimal control problems encountered so far (i.e., up through state 
inequality constraints). We further plan to document the methodology through the completion of 
one paper (which we are now revising in response to reviewers) on the application of the method 
to launch vehicle trajectory analysis, two technical notes on control and state inequality constraints, 
one paper on the general code, and a user’s manual for the code. We have received several calls 
from parties interested in application of the methodology in industry and, although there is nothing 
concrete established as yet, hope to somehow transfer the technology to an industry application in 
the future. 
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Figure 5. lst-order Formulation (Spherical Earth & Back-pressure) 
Checking Using The Hamiltonian H 
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Figure 6. lst-order Formulation (Spherical Earth, Back-pressure & 
Aerodynamic Effects) Checking Using H 
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Figure 7. Angle -of-attack Open-loop Solution Comparison For 
Spherical Earth & Back-pressure Effects 



Figure 8. Flight Path Angle Open-loop Solution Comparison 
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Figure 9. Lamda-v Open-loop Solution Comparison 



Figure 10. Lamda-u Open-loop Solution Comparison 
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Figure 11. Lamda-r Open-loop solution Comparison 
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Appendix A: State Inequality Constraints 

Consider a system defined by a set of n states x and a set of m controls u. Let the 
system be governed by a set of state equations of the form x = f(x,u,t). The class of 
problems to be considered is limited to the case where x is continuous, but there may be 
discontinuities in the costates A. These discontinuities may be the result of state constraints 
present in the problem. Elements of the performance index, J 0 , may be denoted by an 
integrand L(x, u , t) and discrete functions of the states and time 4>[x(t), t] at the initial and 
final times <o and tf. In addition, any constraints imposed on the states and time at the 
initial and final times may be placed in sets of functions rp[x(t), <]. These constraints may 
be adjoined to the performance index by discrete Lagrange multipliers v defined at t 0 and 
tf. Similarly, the state equations may be adjoined to the performance index with a set of 
Lagrange multiplier functions A(t) which will be referred to as costates. 

Now, suppose that there is a scalar constraint on the states and time defined by 
S(x,t ) < 0. The constraint is said to be of p^-order if the p^ total time derivative of 
S is the first to contain the control u explicitly. The first attempt to solve problems with 
state inequality constraints was to use the necessary conditions presented in [12]. These 
necessary conditions lead to successful and accurate solution strategies for states that only 
touch (i.e., do not ride) the constraint boundary. As is derived in [12], for constraints of 
odd order greater than one, the solution can at most only touch the constraint boundary. 
However, for cases where the states ride the constraint boundaries for a nonzero length of 
time, the algebraic equations developed by the weak form are singular. Private discussions 
with Jason Speyer and Dan Moerder indicate that the cause is related to a reduced- 
dimensional manifold; however, we have not been able to develop a nonsingular weak form 
as of now. 

Fortunately, the necessary conditions presented in [8] are accurate for first and second 
order constraints where the solution often rides the constraint boundary. Thus a weak 
formulation is also developed using these necessary conditions for constraints where p = 1 
or p = 2. Therefore, below are presented two very similar weak formulations which are 
accurate for up to a third order constraint and odd-ordered constraints beyond that. Most 
practical applications will be third-order or less. 

General Development 

The weak formulation is now derived for touch-point cases. Without loss of generality, 
assume that there is only one touch-point over the time interval of interest. In this case, 
the state constraint is nothing more than an interior boundary point which creates a jump 
in the costate. 

The performance index J 0 now takes the form: 


f tx f tj 

Jo= [L(x,u,t) + \ T (f — x)] dt+ [L(x,u,t.) + A r (/-i)] dt + viS\ tl +$|!' (1) 
** 1 0 J i 2 
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where $ - <j>[x(t),t] + i/ T rp[x(t),t\. The constraints to be adjoined to J 0 above to trans- 
form the strong boundary conditions to weak boundary conditions are that the states be 
continuous at the initial and final times. Introducing 

x\ to = lim x(t ) and x\ lf = lim x(<) (2) 

i ~ ¥i o t—*tj 

and 


xq = x\ to - x(t 0 ) and xj = x\ t/ = x(t f ) (3) 

one can weakly enforce continuity by adjoining a T (x - x)\\ f Q to J 0 where a is a set of 
discrete unknown Lagrange multipliers defined only at to and tj. The new performance 
index is 


J = Jq + a T (x — 



( 4 ) 


To derive the weak principle, it is necessary to take the first variation of J and set it 
equal to zero. For notations! convenience, the following variables are introduced. 


__ 5$ 

dx 


to 


and A f = 

i dx 


( 5 ) 


Also, as is shown in [4], the Lagrange multiplier a can be chosen so that 8a = 8X. The final 
foim of the weak principle is obtained after integrating by parts so that no derivatives of 
the states or costates appear. After defining the Hamiltonian H = L + A T f and denoting 
the variations of the variables at the initial, touch-point, and final times with subscripts 
0, 1, and / respectively, then the resulting equation is 


/“ I -ti T A + SX T f + SX T x + Sx T \ ( X 
Jt o [ \OX J \OX J 


+ Su T 

+ 8x t 






(It -J - 


f 

Jt 1 


-6i T A + S\ T f + 8\ t x 


+Su r\(m‘ + m‘x 

du J \du J 


dt 


T 

+ 8v'fS\ u +6v T ip\ t ,+6x'[ v\ + SxJX f -Sx^Xo-SXjxf + SXjxo 


+ St} 




dS 

dt 


+ 8t f H (t. f') = 0 


( 6 ) 
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This is the governing equation for the weak Hamiltonian method for problems with touch- 
point state inequality constraints. It is easily shown by integrating the 6x and <5A terms 
by parts in Eq. ( 6 ) that all the Euler-Lagrange equations are the same as in [ 12 ] and that 
all boundary conditions are now of the natural type. 

One simplification may be made to Eq. ( 6 ). If the control is continuous across t\ (as is 
often the case), then it is possible to simplify the 8t\ equation since then /(t]~) = f(tf) = 
f(t l) and £(<f ) = L(t{) = L{ti). From the necessary conditions that are found in [ 12 ] or 
from the ones that could be found from Eq. ( 6 ), it is seen that 

A T (*r)-A T (*+) = i/i^ (7) 


Now, rewriting the coefficient of 8t\ as 


( 3 ) 


we see that the condition for continuity of the Hamiltonian reduces to the condition that 
the first total time derivative of the constraint be zero at t\ if the control is continuous. 

For cases where there is a boundary arc (he., the solution rides the constraint boundary 
for a nonzero length of time), then the weak formulation must be modified. For simplicity 
and without loss of generality, consider the case where the solution has an unconstrained 
arc followed by a constrained arc and then another unconstrained arc. Introducing a new 
Lagrange multiplier function 77 to adjoin the derivative of the constraint S to the 
performance index, then Jo becomes 




dS 

dt 


[A r (ir)-A T ((+)]/(( 1 ) + *,|? 


dS 


dS 


dS 




Jo = / [L(x,u,t) + \ T (f - x)] dt + [ 

Jt 0 Ju 

f tj 

+ / [L(x,u,t) + A T (f - .i-)] dt + + <k|^ 

J t o 


L(x, u, t ) + A T (f — x) + r/ 


dPS' 

dtP 


dt 


(9) 


where N is a column matrix defined as 


N 


T 


S & 


d p * 


dt 


dtP 


4f] 


Again, we define 


( 10 ) 


J = Jo + cd\x-x)\\[ 


( 11 ) 
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Analogous steps to those described above lead to a weak formulation for state constraint 
problems which ride the constraint boundary. These details are presented in [6]. 


Example 

This example is taken from section 3.11 of [8]. The problem is to minimize 



The state equations are 


ii = u 

x 2 — X\ 


(13) 


The state inequality constraint S(x , t) = x 2 - 0 < 0 is to be imposed. For certain values of 
0, the solution only touches the boundary, whereas for other values of 0 the solution rides 
the boundary. 

The algebraic equations were solved using a Newton-Raphson method and a FOR- 
TRAN code written on a SUN 3/260. The sparse, linearized equations are solved using 
subroutine MA28 from the Harwell subroutine library [10]. This subroutine takes advan- 
tage of sparsity which leads to great computational savings. 

The state x 2 is shown in Fig. 12 for the single touch-point case. Results for 2, 4, and 
8 elements on either side of the touch-point (denoted by 2:2, etc.) are compared to the 
exact solution. Note that even the 2:2 element case lies essentially on the exact solution. In 
Fig. 13, the state x 2 is shown for an example case where the state rides the boundary. Here, 
there are three time intervals and the number of elements in each interval is denoted by 
2:2:2 etc. Again we see that the 2:2:2 case has essentially converged on the exact solution. 
One drawback of the weak formulation is that two separate codes had to be written to 
solve this problem. Also, one must determine in advance if the solution will ride or just 
touch the constraint. However, with the general code described in Appendix B, these are 
simple and quick things to do. 
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Figure 13: Displacement vs. Time for a boundary arc case of 0=0. 1 
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Appendix B: General Code Usage 

As an example of how the general code is used, consider the following model of a 
single-stage, four-state rocket. The four states arc m (mass), h (height), V (velocity), and 
7 (flight-path angle). The control u will be the angle-of-attack. Letting T vac be the thrust 
in a vacuum, D be the drag, L be the lift, g be the acceleration due to gravity, I 9p be the 
specific impulse, p be the earth’s gravitational constant, and R e be the radius of the earth, 
then the following equations of motion may be used. 


— T 

■*- V 


m = 


gL 


up 


h = V sin 7 

V = r ~ D _ tL sin 2. 
m (7?. e + h ) 2 

(T + L)u ( V 

1 mV \R e + h 


) 

( Re + hyv ) 


cos 7 


(14) 


For simplicity in this example, the atmospheric pressure has been neglected and the drag 
and lift coefficients have been made constants. Note that this is not necessary in general. 
Thus, 


T = T vac = S155S00N 
p = 1.225 exp( — ft/6700) 

<1 = \pV 2 

D = qS(C Do + Cj<r a U 2 ) 

L = qSC La ( 15 ) 

5 = 33.468 m 2 
C Do = 0.02 
Cno = 6.0 
C La = 5.9S 

The physical constants used in the abov r e model are /t = 3.9906 x 10 14 m 3 s -2 , R e = 6378000 
m, g = 9.81 ms -2 , and I sp = 263.4 s. The performance index is the final mass. The known 
initial conditions are m(0) = 520000 kg, /t(0) = 1800 m, F(0) = 300 m/s, and 7 ( 0 ) = 1.5 
rad. The final conditions arc h(tf ) = 50000 m, V(tf) = 4000 m/s, and 7 (tf) = 0.0 rad. 

Below the input file used to solve this problem is given. The user is required to supply 
the number of states NS, the number of control constraints NP (zero in this example), 
the number of controls M, and the number of constraints on the states at the final time 
Q. The next series of lines from TVAC to F[4] define the system equations as given in 
Eqs. (14) and (15) above. (The lines from TVAC to DRAG are not required but are used 
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to simplify the actual expressions for F[l] thru F[4].) After the equations are formed, the 
user supplies the performance index L and PHI. Then the Q constraints are given in PSI 
and the initial conditions are given in IC. Next the user supplies the final time TF and a 
guess at the value of the final time TFGUES. Since the final time is unknown, TF is set 
to zero and the user gives a guess at the final time. Also, guesses for the states at the 
midpoint of the trajectory and the final point are given in XGUES. These guesses may 
be very crude and can even be zero for many problems. Since the final value of three of 
the states were known for this problem, crude guesses were easily and obviously obtained. 
Finally, the number of elements to be run is given in NE. 

Regardless of the value of NE, the code automatically starts with the two element 
case and uses the continuation method of [11] and the Newton- Raphson method to solve 
the problem. The code then interpolates the solution to this case and runs a four element 
case using only the Newton-Raphson method. The code continues in this manner until 
NE is met. If the Newton-Raphson fails to converge for the four or higher element case 
(which is rare) then the program will start that case over and try the continuation method 
to solve the four element case. 

The output of the example is given after the input file and consists of the solutions 
for the states, costates, controls, and Hamiltonian for 2, 4, and 8 elements. At the top of 
each page is the total elapsed computer time from the start of the program. On the two 
element case sheets is 15.74 secs. This is the time the code took to run the continuation 
method and the Newton-Raphson method for this case. This is a rather small number 
given the complexity of the problem and the fact that an accurate second-order Runge- 
Ivutta method was used to solve the problem. The time at the top of the four element 
case is 18.84 which tells us that only IS. 84 - 15.74 = 3.1 seconds was required to run the 
four-element case given the solution to the two element case. Finally, the desired eight 
element case solution was obtained in a total of 23.46 secs and only 4.62 secs from the 
four element case. Note that this time includes the extraction of nodal values and the 
production of the data files. This is a nonnegligible part of the total time. 

In summary, a complicated rocket trajectory optimization problem which originally 
took several weeks to program and solve is now solved in about 10 or 15 minutes. The 
simple input file is typed in a few minutes and a few minutes are required by MACSYMA 
to create the FORTRAN subroutines. After that, the program runs in a matter of seconds. 
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NS : 4 ; 

NP : 0 ; 

M : 1 ; 

Q : 3 ; 

TVAC: 8155800 . 0; 

ISP : 2 63.4; 

GRAV: 9 . 81; 

MU : 3 . 9906E14 ; 

RE: 6378000.0; 

H (X) : = RE+X (2 ) ; 

RHOSEA: 1.225; 

S : 3 3 . 4 6 8 ; 

CNA : 6 . 0 ; 

CAT : 0 . 0 2 ; 

CLA : 5 . 98 ; 

RHO (X) : =RHOSEA*EXP ( -X ( 2 ) / 67 00 . 0 ) ; 

DP (X) :=0 . 5*RHO (X) *X(3) A 2; 

LIFT (X) : =DP (X) *S*CLA; 

DRAG (X, U) :=DP (X) *S* (CAT+CNA*U (1 ) A 2 ) ; 

F [ 1 ] : -TVAC/ (GRAV* ISP) ; 

F [ 2 ] :X (3) *SIN (X (4 ) ) ; 

F [3] : (TVAC-DRAG (X,U) ) /X (1) - MU*SIN (X (4 ) ) /H (X) A 2; 

F [4] : (TVAC+LIFT (X) ) *U (1 ) / (X (1 ) *X (3) ) + (X (3) /H(X) -MU/ (X (3) *H (X) **2) ) *COS (X ( 4 ) ) 

L : 0 . 0 ; 

PHI :X (1) ; 

PSIfl] :X(2) - 50000.0; 

PSI [2] : X (3) -4000.0; 

PSI [3] : X (4 ) ; 

IC[1] : 520000 . 0 ; 

IC[2] : 1 8 0 0 . 0 ; 

IC [ 3] :300 .0; 

IC [ 4 ] : 1 . 5 ; 

TF : 0 . 0 ; 

TFGUES : 100.0; 

XGUES [1, 1] :260000.0; 

XGUES [1,2] : 10 00.0; 

XGUES [2,1] : 25000. 0; 

XGUES [2,2] : 5 0 0 0 0 . 0 ; 

XGUES [3,1] : 2000 . 0; 

XGUES [3,2] : 4 0 0 0 . 0 ; 

XGUES [ 4 , 1 ] : 0 . 5 ; 

XGUES [ 4 , 2 ] : . 0 ; 

NE : 8 ; 
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NODAL VALUES FOR THE STATES 


NUMBER OF ELEMENTS = 2 TOTAL ELAPSED TIME = 


XI 

0 .52000E+06 
0 .30215E+06 
0 . 84297E+05 


X2 

0 . 18000E+04 
0 . 37 957E+05 
0 . 50000E+05 


X3 

0.30000E+03 
0 . 11357E+04 
0 . 40000E+04 


X4 

0 . 15000E+01 
0 . 13600E+00 
- . 27756E-16 


15.74 


TIME 

0 . 00000E+00 
0.69021E+02 
0 . 13804E+03 
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NODAL VALUES FOR THE STATES 


NUMBER OF ELEMENTS = 4 TOTAL ELAPSED TIME = 


XI 

0 . 52000E+06 
0 . 41393E+06 
0.30786E+06 
0 . 20179E+06 
0 . 95716E+05 


X2 

0 . 18000E+04 
0 . 14334E+05 
0 . 26293E+05 
0 . 40027E+05 
0 .50000E+05 


X3 

0 .30000E+03 
0 . 578 67E+03 
0 . 12035E+04 
0.21967E+04 
0 . 40000E+04 


X4 

0 . 15000E+01 
0 . 52797E+00 
0 .29367E+00 
0 . 19185E+00 
- . 83267E-16 


18 .84 


TIME 

0 .00000E+00 
0.33606E+02 
0 . 67212E+02 
0 . 10082E+03 
0 . 134 42E+03 
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NUMBER 


XI 

0 . 52000E+06 
0 . 4 6737E+06 
0 . 41474E+06 
0.36211E+06 
0.30948E+06 
0.25685E+06 
0.20422E+06 
0.15159E+06 
0 . 98964E+05 


NODAL VALUES FOR THE STATES 


OF ELEMENTS = 8 TOTAL ELAPSED TIME 


X2 

0 .18000E+04 
0 . 71167E+04 
0.12224E+05 
0 . 17557E+05 
0.23577E+05 
0 . 30505E+05 
0.38270E+05 
. 45939E+05 
. 50000E+05 


X3 

0 . 30000E+03 
0 . 39732E+03 
0 . 604 61E+03 
0.88158E+03 
0.12259E+04 
0.16554E+04 
0.22047E+04 
0.29391E+04 
0 . 40000E+04 


X4 

0 . 15000E+01 
0 . 808 61E+00 
0 . 50716E+00 
0 . 38264E+00 
0.31682E+00 
0 . 26830E+00 
0.21910E+00 
0 . 14049E+00 
-.24980E-15 


23.46 


TIME 

0 . OOOOOE+OO 
0 . 16674E+02 
0 . 33349E+02 
0 . 50023E+02 
0 . 66697E+02 
0 . 83371E+02 
0 . 10005E+03 
0.11672E+03 
0 . 13339E+03 
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ALL VALUES FOR CONTROL AND HAMILTONIAN 


NUMBER OF ELEMENTS = 2 TOTAL ELAPSED TIME = 


U1 

- . 92485E+00 
- . 27303E+00 
0.42760E+00 
0 . 74199E-01 
0.33028E-01 


U2 

0 .00000E+00 
0 .00000E+00 
0 . 00000E+00 
0 . 00000E+00 
0 . 00000E+00 


U3 

0 . 00000E+00 
0 . 00000E+00 
0 . 00000E+00 
0 . 00000E+00 
0 . 00000E+00 


HAMIL 

- . 4 6274E+03 
- . 82633E+03 
- . 71640E+03 
98340E+03 
- . 10640E-08 


15 .74 


TIME 

0 . 00000E+00 
0 .34510E+02 
0 . 69021E+02 
0 .10353E+03 
0 . 13804E+03 
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ALL VALUES FOR CONTROL AND HAMILTONIAN 


NUMBER OF ELEMENTS = 4 TOTAL ELAPSED TIME = 


U1 

- . 61026E+00 
- . 23120E+00 
- . 12248E-01 
0.73313E-01 
0 . 14683E+00 
0 . 91624E-01 
0 . 56181E-01 
14828E+00 
- . 32205E+00 


U2 

0 .00000E+00 
0 . 00000E+00 
0 . 00000E+00 
0 . OOOOOE+OO 
0 . 00000E+00 
0 .00000E+00 
0 . 00000E+00 
0. 00000E+00 
0 . 00000E+00 


U3 

0 .00000E+00 
0 . OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 
0. 00000E+00 
0. OOOOOE+OO 
0 .00000E+00 
0 . OOOOOE+OO 
0 . OOOOOE+OO 


HAMIL 

- . 17105E+03 
- . 288 90E+03 
- . 23618E+03 
- . 24606E+03 
- . 23657E+03 
-.26451E+03 
- . 21651E+03 
- . 38816E+03 
0 . 24727E-11 


18.84 


TIME 

0. OOOOOE+OO 
0 .16803E+02 
0 .33606E+02 
0 . 50409E+02 
0.67212E+02 
0 . 84015E+02 
0.10082E+03 
0 . 11762E+03 
0 . 13442E+03 
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ALL VALUES FOR CONTROL AND HAMILTONIAN 


NUMBER OF ELEMENTS = 8 TOTAL ELAPSED TIME = 


U1 

- . 55186E+00 
- . 32779E+00 
- . 16595E+00 
- . 38584E-01 
0 . 44739E-01 
0 . 79906E-01 
0 . 10973E+00 
0 . 11215E+00 
0 . 12021E+00 
0.10770E+00 
0 . 10186E+00 
0 . 67127E-01 
0 .29691E-01 
66484E-01 
- . 18969E+00 
- . 28902E+00 
- . 37039E+00 


U2 

0 .00000E+00 
0 .00000E+00 
0 .OOOOOE+OO 
0.00000E+00 
0 .OOOOOE+OO 
0 .00000E+00 
0 . OOOOOE+OO 
0. 00000E+00 
0 . 00000E+00 
0 .OOOOOE+OO 
0 .00000E+00 
0 .OOOOOE+OO 
0. OOOOOE+OO 
0 .00000E+00 
0. OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 


U3 

0 .OOOOOE+OO 
0. OOOOOE+OO 
0 .OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 
0 .OOOOOE+OO 
0. OOOOOE+OO 
0 .OOOOOE+OO 
0 . OOOOOE+OO 
0 .OOOOOE+OO 
0 .OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 
0. OOOOOE+OO 


HAMIL 

- . 51935E+02 
- . 95367E+02 
- . 62856E+02 
- . 764 63E+02 
- . 69769E+02 
- . 71521E+02 
- . 70243E+02 
- . 72219E+02 
- . 69462E+02 
- . 74437E+02 
- . 67925E+02 
- . 79036E+02 
- . 64306E+02 
- . 92023E+02 
- . 52356E+02 
- . 14072E+03 
0 . 16485E-10 


i 


23.46 


TIME 

0 .OOOOOE+OO 
0 .83371E+01 
0 . 16674E+02 
0 . 250 11E+02 
0.33349E+02 
0 .41686E+02 
0 .50023E+02 
0 . 58360E+02 
0 . 66697E+02 
0 . 75034E+02 
0 .83371E+02 
0 .91709E+02 
0 .10005E+03 
0 .10838E+03 
0 .11672E+03 
0.12506E+03 
0.13339E+03 
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NODAL VALUES FOR THE COSTATES 


NUMBER OF ELEMENTS = 2 TOTAL ELAPSED TIME = 


LI 

0 . 38777E+00 
0.50266E+00 
0 . 10000E+01 


L2 

0 . 238 95E+00 
0 .13755E+00 
0 . 12191E+00 


L3 

0 .33772E+02 
0 .33298E+02 
0.32653E+02 


L4 

95589E+04 

0.20419E+04 

0.10502E+04 


15.74 


TIME 

0 . 00000E+00 
0 . 69021E+02 
0 . 13804E+03 
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NODAL VALUES FOR THE COSTATES 


NUMBER OF ELEMENTS = 4 TOTAL ELAPSED TIME = 


LI 

0 .24671E+00 
0 .30216E+00 
0 . 38809E+00 
0 . 54729E+00 
0 . 10000E+01 


L2 

0 . 31034E+00 
0 . 28438E+00 
0 . 26974E+00 
0 .25414E+00 
0 . 22100E+00 


L3 

0 . 41713E+02 
0 . 42814E+02 
0.38756E+02 
0 .36560E+02 
0 .36298E+02 


L4 

- . 77904E+04 
- . 22655E+03 
0.41327E+04 
0 . 14102E+04 
- . 11383E+05 


18 .84 


TIME 

0 .00000E+00 
0 . 33606E+02 
0 .67212E+02 
0 . 10082E+03 
0.13442E+03 
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NODAL VALUES FOR THE COSTATES 


NUMBER OF ELEMENTS = 8 TOTAL ELAPSED TIME = 


LI 

0 .19901E+00 
0 .22653E+00 
0.25878E+00 
0.29819E+00 
0 .34844E+00 
0.41599E+00 
0 . 51405E+00 
0.67453E+00 
0 . 10000E+01 


L2 

0 .34293E+00 
0 . 32947E+00 
0 .32159E+00 
0 . 312 94E+00 
0 .30244E+00 
0 .29121E+00 
0 . 28441E+00 
0 . 27907E+00 
0 . 22278E+00 


L3 

0 . 42300E+02 
0 .47525E+02 
0 .44891E+02 
0.42327E+02 
0.40402E+02 
0 . 38903E+02 
0 . 37713E+02 
0 .36970E+02 
0 . 37315E+02 


L4 

— . 714 41E+04 
- . 28336E+04 
0 . 11449E+04 
0.37755E+04 
0.47900E+04 
0 . 39834E+04 
0 . 96373E+03 
- . 4 9726E+04 
- . 13459E+05 


23 .46 


TIME 

0 .00000E+00 
0 .16674E+02 
0 .33349E+02 
0 .50023E+02 
0 . 66697E+02 
0 . 83371E+02 
0 . 100 05E+03 
0 .11672E+03 
0 . 13339E+03 
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